Take-home_Ex1

Introduction

The success of Singapore’s public transport is commonly recognised as 1 of the quintessential indicator of its rapid growth and prosperity since her independence in 1965. It comprises of both public buses as well as mass rapid transport (MRT) trains. While the MRT network has been expanding across the many zones within Singapore, the bus network remains the most option accessible and reliable option amongst the two.

Thus, the utilisation pattern of buses are of key importance to multiple aspects of Singapore’s socio-economic spheres, ranging from areas of work productivity and efficiency, environmental impact, tourism, and even to potential impact to cost of living and real estate development; which underpins the impetus of this exercise to gain a better understanding of the utilisation patterns and characteristics of Singapore’s bus stops.

Objectives

The key objectives of this exercise are:

• To perform geovisualisation and analysis of peakhour busstop passenger trips
• To perform Analysis of Local Indicators of Spatial Association (LISA) for passenger trips

Getting Started

The code chunk below load the following packages:

• tmap: for thematic mapping
• sf for geospatial data handling
• tidyverse for non-spatial data handling.
• knitr for creating html table
• spdep
• sfdep for creating an sf and tidyverse friendly interfaceas well as for introducing new functionality that is not present in spdep
• plotly for interactive plots

pacman::p_load(tmap, sf, tidyverse, knitr, spdep, sfdep, plotly, zoo, Kendall)

Importing and preparing the Base Geospatial and Aspatial Data

Aspatial data

Firstly, the Passenger Volume by Origin Destination Bus Stops data set for the month of October 2023, from LTA DataMall, will be imported by using ‘read_csv()’ of **readr* package.

#\| eval: false 
odbus <- read.csv("data/aspatial/origin_destination_bus_202310.csv")

From below, the values in OROGIN_PT_CODE and DESTINATON_PT_CODE odbus in the tibble data frame are in numeric data type.

glimpse(odbus)

Rows: 5,694,297
Columns: 7
$ YEAR_MONTH          <chr> "2023-10", "2023-10", "2023-10", "2023-10", "2023-…
$ DAY_TYPE            <chr> "WEEKENDS/HOLIDAY", "WEEKDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <int> 16, 16, 14, 14, 17, 17, 17, 7, 14, 14, 10, 20, 20,…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <int> 4168, 4168, 80119, 80119, 44069, 20281, 20281, 190…
$ DESTINATION_PT_CODE <int> 10051, 10051, 90079, 90079, 17229, 20141, 20141, 1…
$ TOTAL_TRIPS         <int> 3, 5, 3, 5, 4, 1, 24, 2, 1, 7, 3, 2, 5, 1, 1, 1, 1…

Using appropriate tidyverse functions to convert these data values into factor data type.

odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE)

From below, the values in OROGIN_PT_CODE and DESTINATON_PT_CODE odbus in the tibble data frame are now in factor data type.

glimpse(odbus)

Rows: 5,694,297
Columns: 7
$ YEAR_MONTH          <chr> "2023-10", "2023-10", "2023-10", "2023-10", "2023-…
$ DAY_TYPE            <chr> "WEEKENDS/HOLIDAY", "WEEKDAY", "WEEKENDS/HOLIDAY",…
$ TIME_PER_HOUR       <int> 16, 16, 14, 14, 17, 17, 17, 7, 14, 14, 10, 20, 20,…
$ PT_TYPE             <chr> "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "BUS", "…
$ ORIGIN_PT_CODE      <fct> 4168, 4168, 80119, 80119, 44069, 20281, 20281, 190…
$ DESTINATION_PT_CODE <fct> 10051, 10051, 90079, 90079, 17229, 20141, 20141, 1…
$ TOTAL_TRIPS         <int> 3, 5, 3, 5, 4, 1, 24, 2, 1, 7, 3, 2, 5, 1, 1, 1, 1…

Geospatial data

Import BusStop layer downloaded from LTA DataMall into RStudio and save it as a sf data frame called busstop.

busstop <- st_read(dsn = "data/geospatial",
                   layer = "BusStop") %>%
  st_transform(crs = 3414)

Reading layer `BusStop' from data source 
  `C:\jayexx\ISSS624\Take-home_Exercises\Take-home_Ex1\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5161 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

The structure of busstop sf tibble data frame is as follows.

glimpse(busstop)

Rows: 5,161
Columns: 4
$ BUS_STOP_N <chr> "22069", "32071", "44331", "96081", "11561", "66191", "2338…
$ BUS_ROOF_N <chr> "B06", "B23", "B01", "B05", "B05", "B03", "B02A", "B02", "B…
$ LOC_DESC   <chr> "OPP CEVA LOGISTICS", "AFT TRACK 13", "BLK 239", "GRACE IND…
$ geometry   <POINT [m]> POINT (13576.31 32883.65), POINT (13228.59 44206.38),…

Objective (1): Geovisualisation and Analysis of Bus Stop Passengers during Peakhours

Extracting the passenger trips for peakhour time intervals

For the purpose of this exercise, extract commuting flows during:

• weekday morning peak and label the output tibble data table as odbus_wkd6_9.

odbus_wkd6_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 &
           TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

• weekday evening peak and label the output tibble data table as odbus_wkd17_20.

odbus_wkd17_20 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 17 &
           TIME_PER_HOUR <= 20) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

• weekend/holiday morning peak and label the output tibble data table as odbus_wke11_14.

odbus_wke11_14 <- odbus %>%
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 11 &
           TIME_PER_HOUR <= 14) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

• weekend/holiday evening peak and label the output tibble data table as odbus_wke16_19.

odbus_wke16_19 <- odbus %>%
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 16 &
           TIME_PER_HOUR <= 19) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

Convert and save the tible data sets in rds format as follows.

write_rds(odbus_wkd6_9, "data/rds/odbus_wkd6_9.rds")
write_rds(odbus_wkd17_20, "data/rds/odbus_wkd17_20.rds")
write_rds(odbus_wke11_14, "data/rds/odbus_wke11_14.rds")
write_rds(odbus_wke16_19, "data/rds/odbus_wke16_19.rds")

Import files back into R as follows.

odbus_wkd6_9 <- read_rds("data/rds/odbus_wkd6_9.rds")
odbus_wkd17_20 <- read_rds("data/rds/odbus_wkd17_20.rds")
odbus_wke11_14 <- read_rds("data/rds/odbus_wke11_14.rds")
odbus_wke16_19 <- read_rds("data/rds/odbus_wke16_19.rds")

Preparing Geospatial data in Hexagonal grid

Creating hexagonal grid in sf format with grid_id as follows.

hxgn_grid = st_make_grid(busstop, c(250, 250), what = "polygons", square = FALSE)

# To sf and add grid ID
hxgn_grid_sf = st_sf(hxgn_grid) %>%
  mutate(grid_id = 1:length(lengths(hxgn_grid)))

Merging imported busstop data with Hexagonal grid as follows.

busstop_hxgn_grid <- st_intersection(busstop, hxgn_grid_sf) %>%
  select(BUS_STOP_N, grid_id) %>%
  st_drop_geometry()

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Checking for BUS_STOP_N with duplicate grid_id, and grid_id with duplicate BUS_STOP_N.

check_duplicate <- busstop_hxgn_grid %>%
  group_by(grid_id) %>%
  summarise(num_BUS_STOP_N = n_distinct(BUS_STOP_N))

check_duplicate1 <- busstop_hxgn_grid %>%
  group_by(BUS_STOP_N) %>%
  summarise(num_grid_id = n_distinct(grid_id))

From the above, - out of 3,131 grid_id, 1,476 (47%) grid_id contain 2 busstops, 1,410 (45%) grid_id contain only 1 busstop, and the remaining 8% containing more than 2 busstops. Only 1 grid_id has the maximum of 5 busstops. - out of 5,145 busstops, only 8 intercept across 2 grid_id.

Hence, for this exercise, to avoid double counting of trips subsequently, the duplicate grid_id for the 8 busstops are removed.

busstop_hxgn_grid <- distinct(busstop_hxgn_grid, BUS_STOP_N, .keep_all = TRUE)

Merging Peakhour Passenger Trips with Geospatial data

Weekday morning peak

Append the merged busstop_hxgn_grid frame onto odbus_wkd6_9 data frame as follows.

wkd6_9_hxgn_grid <- left_join(busstop_hxgn_grid, odbus_wkd6_9,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) %>%
  rename(ORIGIN_BS = BUS_STOP_N,
         ORIGIN_GRID = grid_id) %>%
  group_by(ORIGIN_GRID) %>%
  reframe(
    MAX_TRIPS_BS = ORIGIN_BS[which.max(sum(TRIPS))],
    TOT_TRIPS = sum(TRIPS),
    AVG_TRIPS = mean(TRIPS)
  ) 

As a good practice, check for duplicate records and retain duplicate records as follows.

duplicate <- wkd6_9_hxgn_grid %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

wkd6_9_hxgn_grid <- unique(wkd6_9_hxgn_grid)

Merge resultant data with hexgon grid and filter as follows

wkd6_9_hxgn_grid_sf <- left_join(hxgn_grid_sf, 
                           wkd6_9_hxgn_grid,
                           by = c("grid_id" = "ORIGIN_GRID"))%>%
  select(MAX_TRIPS_BS, TOT_TRIPS, AVG_TRIPS, hxgn_grid)

wkd6_9_hxgn_grid_sf = filter(wkd6_9_hxgn_grid_sf, TOT_TRIPS > 0)

Weekday evening peak

Append the merged busstop_hxgn_grid frame onto odbus_wkd17_20 data frame as follows.

wkd17_20_hxgn_grid <- left_join(busstop_hxgn_grid, odbus_wkd17_20,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) %>%
  rename(ORIGIN_BS = BUS_STOP_N,
         ORIGIN_GRID = grid_id) %>%
  group_by(ORIGIN_GRID) %>%
  reframe(
    MAX_TRIPS_BS = ORIGIN_BS[which.max(sum(TRIPS))],
    TOT_TRIPS = sum(TRIPS),
    AVG_TRIPS = mean(TRIPS)
  )

As a good practice, check for duplicate records and retain duplicate records as follows.

duplicate1 <- wkd17_20_hxgn_grid %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

wkd17_20_hxgn_grid <- unique(wkd17_20_hxgn_grid)

Merge resultant data with hexgon grid and filter as follows

wkd17_20_hxgn_grid_sf <- left_join(hxgn_grid_sf, 
                           wkd17_20_hxgn_grid,
                           by = c("grid_id" = "ORIGIN_GRID"))%>%
  select(MAX_TRIPS_BS, TOT_TRIPS, AVG_TRIPS, hxgn_grid)

wkd17_20_hxgn_grid_sf = filter(wkd17_20_hxgn_grid_sf, TOT_TRIPS > 0)

Weekend/holiday morning peak

Append the merged busstop_hxgn_grid frame onto odbus_wke11_14 data frame as follows.

wke11_14_hxgn_grid <- left_join(busstop_hxgn_grid, odbus_wke11_14,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) %>%
  rename(ORIGIN_BS = BUS_STOP_N,
         ORIGIN_GRID = grid_id) %>%
  group_by(ORIGIN_GRID) %>%
  reframe(
    MAX_TRIPS_BS = ORIGIN_BS[which.max(sum(TRIPS))],
    TOT_TRIPS = sum(TRIPS),
    AVG_TRIPS = mean(TRIPS)
  )

As a good practice, check for duplicate records and retain duplicate records as follows.

duplicate2 <- wke11_14_hxgn_grid %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

wke11_14_hxgn_grid <- unique(wke11_14_hxgn_grid)

Merge resultant data with hexgon grid and filter as follows

wke11_14_hxgn_grid_sf <- left_join(hxgn_grid_sf, 
                           wke11_14_hxgn_grid,
                           by = c("grid_id" = "ORIGIN_GRID"))%>%
  select(MAX_TRIPS_BS, TOT_TRIPS, AVG_TRIPS, hxgn_grid)

wke11_14_hxgn_grid_sf = filter(wke11_14_hxgn_grid_sf, TOT_TRIPS > 0)

Weekend/holiday evening peak

Append the merged busstop_hxgn_grid frame onto odbus_wke16_19 data frame as follows.

wke16_19_hxgn_grid <- left_join(busstop_hxgn_grid, odbus_wke16_19,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) %>%
  rename(ORIGIN_BS = BUS_STOP_N,
         ORIGIN_GRID = grid_id) %>%
  group_by(ORIGIN_GRID) %>%
  reframe(
    MAX_TRIPS_BS = ORIGIN_BS[which.max(sum(TRIPS))],
    TOT_TRIPS = sum(TRIPS),
    AVG_TRIPS = mean(TRIPS)
  )

As a good practice, check for duplicate records and retain duplicate records as follows.

duplicate3 <- wke16_19_hxgn_grid %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

wke16_19_hxgn_grid <- unique(wke16_19_hxgn_grid)

Merge resultant data with hexgon grid and filter as follows

wke16_19_hxgn_grid_sf <- left_join(hxgn_grid_sf, 
                           wke16_19_hxgn_grid,
                           by = c("grid_id" = "ORIGIN_GRID"))%>%
  select(MAX_TRIPS_BS, TOT_TRIPS, AVG_TRIPS, hxgn_grid)

wke16_19_hxgn_grid_sf = filter(wke16_19_hxgn_grid_sf, TOT_TRIPS > 0)

Geovisualisation of Busstop Peakhour periods

Plot the map with hexagonal grid as follows for:

• Weekday morning peak

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(wkd6_9_hxgn_grid_sf) +
  tm_fill(
    col = c("TOT_TRIPS", "AVG_TRIPS"),
    palette = "Reds",
    style = "quantile",
    title = c("Total Passenger Trips during Weekday Morning Peak", "Average Passenger Trips during Weekday Morning Peak"),
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Total Trips: " = "TOT_TRIPS",
      "Average Trips: " = "AVG_TRIPS",
      "Most Popular Bus Stop: " = "MAX_TRIPS_BS"
      
    ),
    popup.format = list(
      TOT_TRIPS = list(format = "f", digits = 0),
      AVG_TRIPS = list(format = "f", digits = 0),
      MAX_TRIPS_BS = list(format = "f", digits = 0)
    )
  ) +
  tm_layout(main.title = "Bus Stop Passenger Trips during Weekday Morning Peak",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5)

• Weekday evening peak

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(wkd17_20_hxgn_grid_sf) +
  tm_fill(
    col = c("TOT_TRIPS", "AVG_TRIPS"),
    palette = "Reds",
    style = "quantile",
    title = c("Total Passenger Trips during Weekday Evening Peak", "Average Passenger Trips during Weekday Evening Peak"),
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Total Trips: " = "TOT_TRIPS",
      "Average Trips: " = "AVG_TRIPS",
      "Most Popular Bus Stop: " = "MAX_TRIPS_BS"
      
    ),
    popup.format = list(
      TOT_TRIPS = list(format = "f", digits = 0),
      AVG_TRIPS = list(format = "f", digits = 0),
      MAX_TRIPS_BS = list(format = "f", digits = 0)
    )
  ) +
  tm_layout(main.title = "Bus Stop Passenger Trips during Weekday Evening Peak",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5)

• Weekend/holiday morning peak

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(wke11_14_hxgn_grid_sf) +
  tm_fill(
    col = c("TOT_TRIPS", "AVG_TRIPS"),
    palette = "Reds",
    style = "quantile",
    title = c("Total Passenger Trips during Weekend Morning Peak", "Average Passenger Trips during Weekend Morning Peak"),
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Total Trips: " = "TOT_TRIPS",
      "Average Trips: " = "AVG_TRIPS",
      "Most Popular Bus Stop: " = "MAX_TRIPS_BS"
      
    ),
    popup.format = list(
      TOT_TRIPS = list(format = "f", digits = 0),
      AVG_TRIPS = list(format = "f", digits = 0),
      MAX_TRIPS_BS = list(format = "f", digits = 0)
    )
  ) +
  tm_layout(main.title = "Bus Stop Passenger Trips during Weekend Morning Peak",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5)

• Weekend/holiday evening peak

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(wke16_19_hxgn_grid_sf) +
  tm_fill(
    col = c("TOT_TRIPS", "AVG_TRIPS"),
    palette = "Reds",
    style = "quantile",
    title = c("Total Passenger Trips during Weekend Evening Peak", "Average Passenger Trips during Weekend Evening Peak"),
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Total Trips: " = "TOT_TRIPS",
      "Average Trips: " = "AVG_TRIPS",
      "Most Popular Bus Stop: " = "MAX_TRIPS_BS"
      
    ),
    popup.format = list(
      TOT_TRIPS = list(format = "f", digits = 0),
      AVG_TRIPS = list(format = "f", digits = 0),
      MAX_TRIPS_BS = list(format = "f", digits = 0)
    )
  ) +
  tm_layout(main.title = "Bus Stop Passenger Trips during Weekend Evening Peak",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5)

Spatial patterns revealed by the geovisualisation

From the Geovisualisation of the 4 peak periods, it is noted that the number of trips during weekdays is higher than weekends/holidays; while within the same type of days, evening is higher than morning. Amongst the bus stops, the highest number of trips across the 4 peak periods is from the bus interchanges; while by region, the southwestern region of Singapore seems to have the least number of trips across the 4 peak periods.

Analysis of Local Indicators of Spatial Association (LISA)

Merging Passenger Trips with Geospatial data

Append the merged busstop_hxgn_grid frame onto the complete odbus data frame as follows.

odbus_hxgn_grid <- left_join(busstop_hxgn_grid, odbus,
            by = c("BUS_STOP_N" = "ORIGIN_PT_CODE")) %>%
  rename(ORIGIN_BS = BUS_STOP_N,
         ORIGIN_GRID = grid_id) %>%
  group_by(ORIGIN_GRID) %>%
  reframe(
    MAX_TRIPS_BS = ORIGIN_BS[which.max(sum(TOTAL_TRIPS))],
    TOT_TRIPS = sum(TOTAL_TRIPS),
    AVG_TRIPS = mean(TOTAL_TRIPS)
  ) 

As a good practice, check for duplicate records and retain duplicate records as follows.

duplicate4 <- odbus_hxgn_grid %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

odbus_hxgn_grid <- unique(odbus_hxgn_grid)

Merge resultant data with hexgon grid and filter as follows

odbus_hxgn_grid_sf <- left_join(hxgn_grid_sf, 
                           odbus_hxgn_grid,
                           by = c("grid_id" = "ORIGIN_GRID"))%>%
  select(MAX_TRIPS_BS, TOT_TRIPS, AVG_TRIPS, hxgn_grid)

odbus_hxgn_grid_sf = filter(odbus_hxgn_grid_sf, TOT_TRIPS > 0)

Deriving adaptive distance based-weights

derive an adaptive spatial weights by using the code chunk below.

wm_ad <- odbus_hxgn_grid_sf %>% 
  mutate(nb = st_knn(hxgn_grid,
                     k=8),
         wt = st_weights(nb),
               .before = 1)

! Polygon provided. Using point on surface.

Computing local Moran’s I

Compute Local Moran’s I of GDPPC at county level by using local_moran() of sfdep package as follows.

lisa <- wm_ad %>% 
  mutate(local_moran = local_moran(
    TOT_TRIPS, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)

Visualising local Moran’s I and p-value

For effective comparison, it will be better for us to plot both maps next to each other as shown below.

tmap_mode("view")

tmap mode set to interactive viewing

map1 <- tm_shape(lisa) +
  tm_fill("ii") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "local Moran's I of Passenger Trips",
            main.title.size = 0.8)

map2 <- tm_shape(lisa) +
  tm_fill("p_ii_sim",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig")) + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of local Moran's I",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

Variable(s) "ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Visualising LISA map

LISA map is a categorical map showing outliers and clusters. There are two types of outliers namely: High-Low and Low-High outliers. Likewise, there are two type of clusters namely: High-High and Low-Low cluaters. In fact, LISA map is an interpreted map by combining local Moran’s I of geographical areas and their respective p-values.

In lisa sf data.frame, we can find three fields contain the LISA categories. They are mean, median and pysal. In general, classification in mean will be used as shown in the code chunk below.

lisa_sig <- lisa  %>%
  filter(p_ii_sim < 0.05)
tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(lisa) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(lisa_sig) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).

Prepare a choropleth map showing the distribution of passenger trips at planning sub-zone level as follows

tm_shape(origintrip_SZ)+ tm_fill(“TOT_TRIPS”, style = “quantile”, palette = “Blues”, title = “Passenger trips”) + tm_layout(main.title = “Passenger trips generated at planning sub-zone level”, main.title.position = “center”, main.title.size = 1.2, legend.height = 0.45, legend.width = 0.35, frame = TRUE) + tm_borders(alpha = 0.5) + tm_compass(type=“8star”, size = 2) + tm_scale_bar() + tm_grid(alpha =0.2) + tm_credits(“Source: Planning Sub-zone boundary from URAand Passenger trips data from LTA”, position = c(“left”, “bottom”))

Creating interactive map

tmap_mode(“view”) tmap_options(check.and.fix = TRUE) tm_shape(origintrip_SZ)+ tm_fill(“TOT_TRIPS”, style = “quantile”, palette = “Blues”, title = “Passenger trips”) + tm_layout(main.title = “Passenger trips generated at planning sub-zone level”, main.title.position = “center”, main.title.size = 1.2, legend.height = 0.45, legend.width = 0.35, frame = TRUE) + tm_borders(alpha = 0.5) + tm_compass(type=“8star”, size = 2) + tm_scale_bar() + tm_grid(alpha =0.2) + tm_credits(“Source: Planning Sub-zone boundary from URAand Passenger trips data from LTA”, position = c(“left”, “bottom”))